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This paper reports results from numerical simulations of the gravitational radiation emitted from 
non-rotating compact objects (both neutron stars and Schwarzschild black holes) as a result of the 
accretion of matter. We adopt a hybrid procedure in which we evolve numerically, and assuming 
axisymmetry, the linearized equations describing metric and fluid perturbations coupled to a fully 
nonlinear hydrodynamics code that calculates the motion of the accreting matter. The initial 
matter distribution, which is initially at rest, is shaped in the form of extended quadrupolar shells 
of either dust or obeying a perfect fluid equation of state. Self-gravity of the accreting layers 
of fluid is neglected, as well as radiation reaction effects. We use this idealized setup in order 
to understand the qualitative features appearing in the energy spectrum of the gravitational wave 
emission from compact stars or black holes, subject to accretion processes involving extended objects. 
A comparison for the case of point-like particles falling radially onto black holes is also provided. 
Our results show that, when the central object is a black hole, the spectrum is far from having 
only one clear, monochromatic peak at the frequency of the fundamental quasi-normal mode. On 
the contrary, it shows a complex pattern, with distinctive interference fringes produced by the 
interaction between the infalling matter and the underlying perturbed spacetime, in close agreement 
with results for point-like particles. Remarkably, most of the energy is emitted at frequencies lower 
than that of the fundamental mode of the black hole. Similar results are obtained for extended 
shells accreting onto neutron stars, but in this case the contribution of the stellar fundamental 
mode stands clearly in the energy spectrum. Our analysis illustrates that the gravitational wave 
signal driven by accretion onto compact objects is influenced more by the details and dynamics of 
the process, and the external distribution of matter, than by the quasi-normal mode structure of 
the central object. The gravitational waveforms from such accretion events appear to be much more 
complex than former simplifled assumptions predicted. 

PACS numbers: 04.30.Db, 04.40.Dg, 95.30.Lz, 98.62.Mw 



I. INTRODUCTION 

Most of the different astrophysical scenarios suggested 
as potential sources of gravitational radiation have in 
common the presence of compact objects, such as neu- 
tron stars, strange stars or black holes. The coales- 
cence of a binary system formed by two black holes, 
two neutron stars or one black hole and one neutron 
star is the main target of the ground based interferom- 
eters [Laser Interferometer Gravitational Wave Observer 
(LIGO), VIRGO), but other possibilities such as galactic 
supernovae are also worth exploring, especially in antic- 
ipation of the capabilities of future detectors. Isolated 
black holes, in particular, are reckoned to be character- 
ized by a unique emission pattern known as quasi-normal 
mode (QNM) ringing-rapidly damped sinusoidal modes. 
This signal has been studied extensively using perturba- 
tion theory and frequency-domain techniques for most 
classes of black hole solutions (see e.g. P, Q and ref- 
erences therein). The detection of such QNM signals 
depends strongly on the luminosity of the source or, in 
other words, on how strongly the black hole is excited, 
but also on the knowledge of the power spectrum. It is 
also well known that a relativistic star has a very rich 



non-radial oscillation spectrum, and it can emit gravi- 
tational waves through QNM ringing (see e.g. [J and 
references therein). Therefore, another plausible source 
of gravitational radiation involves the rapid accretion of 
large clumps of matter onto a compact object (either a 
neutron star or a black hole), a recurrent and ubiquitous 
phenomenon in relativistic astrophysics. Accretion is ex- 
pected to happen following the gravitational collapse of 
the core of a massive star, once a neutron star has al- 
ready been formed. Part of the remaining stellar mate- 
rial, which has not been expelled by the shock driving 
the supernova explosion, may fall back onto the neutron 
star, until a critical mass is exceeded and the star col- 
lapses to a black hole. Some more material may in turn 
form a long-lived, centrifugally supported torus or disk if 
the collapsing star had initially some amount of rotation. 

Semi-analytical studies of extended objects (shells or 
blobs of dust) falling isotropically onto a black hole are 
available in the literature |E |1 |l H H . These stud- 
ies shed the first light in understanding the modification 
of the gravitational wave (GW) emission pattern of the 
black hole due to the presence of matter. Collectively, 
these works showed that for a fixed amount of infalling 
mass m, the energy released in gravitational radiation is 
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reduced compared to the value of the point -particle limit 
{E O.Olm^/Af ^, M being the mass of the black hole). 
This reduction is interpreted as due to cancellations of 
the emission from distinct parts of the extended object. 
Such conclusions were later confirmed by Papadopoulos 
and Font 10], who performed numerical simulations of 
the gravitational radiation emitted during the accretion 
process of an extended object onto a black hole. In ^3 
the first-order deviations from the exact black hole geom- 
etry were approximated by the Teukolsky equation [Tll | 
for Schwarzschild black holes, i.e. the inhomogeneous 
Bardeen-Press equation 'l2| , including curvature pertur- 
bations induced by matter sources, whose nonlinear evo- 
lution was integrated using a hydrodynamics code. This 
was the first numerical study in the time-domain of the 
gravitational radiation emitted by extended objects ac- 
creted by black holes, and showed the gradual excitation 
of the black hole fundamental QNM frequency by suffi- 
ciently compact shells. In the thin shell limit, the energy 
asymptotes to a finite value, which is about a third of 
the point-particle limit. 

Correspondingly, linear perturbation studies in the 
time domain of neutron star spacetimes, aimed at ad- 
dressing QNM excitation, have also received consider- 
able attention in the literature (see e.g. and refer- 
ences therein). The oscillation properties are obtained 
from the analysis of non-radial stellar perturbations. For 
non-rotating stars with polytropic equations of state, the 
spectrum naturally splits into an axial (or odd-parity) 
and a polar (or even-parity) part, according to the termi- 
nology used to address the behaviour of the perturbation 
equations under parity transformations. The axial part 
of the spectrum includes only gravitational modes, the 
so called w modes, which a re p urely relativistic, being 
absent in Newtonian gravity |13| . The polar part, on the 
other hand, contains w modes as well as fluid modes. The 
excitation of axial parity modes was analyzed by Anders- 
son and Kokkotas [iJI i by sending pulses of gravitational 
waves to the neutron star and studying its response, and 
by Ferrari and Kokkotas by scattering of point-like 
particles. The excitation of polar modes was first inves- 
tigated by Allen et al. , also as a scattering problem 
using gravitational wave pulses. The same approach was 
later followed by Ruoff and co-workers in T7, us- 
ing both Gaussian pulses and point-particle scattering. 
A framework for constructing initial data sets for per- 
turbations was developed in 19] and applied to study 
neutron star collisions in the close limit approximation 
in Ref. 

Alternately, in a series of papers Seidel and co-workers 
computed the gravitational radiation emanating from 
slightly non-spherical stellar core collapse, in the axial 
plj as well as in the polar case [S^] , and the waveforms as- 
sociated with the formation of neutron stars. The zeroth- 
order solution was a spherical collapsing star, whose dy- 
namics was computed by solving the coupled system 
of Einstein and hydrodynamics equations using the La- 
grangian May- White approach. The GWs were extracted 



using perturbation theory on the spherical background 
within the Gerlach-Sengupta [23| formalism. More re- 
cently, Harada et al. p4|] have reexamined the axial 
part of this problem, using null coordinates (Hernandez- 
Misner) and a gauge invariant and coordinate indepen- 
dent perturbative formalism developed by Marti'n-Garcia 
and Gundlach [2^ |2^ |23 . Within this approach Harada 
et al. [24] have been able to follow the spherical col- 
lapse of both, supermassive stars and neutron stars, until 
a black hole forms, computing the GWs that are emit- 
ted in the process. Full numerical relativity computa- 
tions of QNM ringing in spherical gravitational collapse 
using the characteristic formulation of general relativity 
are also reported in [S^. Within the same formalism, 
the imprints of fluid accretion on the emitted GWs from 
a black hole were studied in .29] . finding the familiar 
damped-oscillatory GW decay, but both decay rate and 
frequencies being modulated by the mass accretion rate. 

In this paper, we analyze the similarities and differ- 
ences in the gravitational wave emission pattern from 
black holes and neutron stars as a result of the radial ac- 
cretion of matter. A detailed realistic modelization of the 
gravitational emission from accretion flows would require 
of three-dimensional (magneto-)hydrodynamical simula- 
tions in general relativity, coupled to radiation trans- 
port and diffusive processes. However, some preliminary 
steps can be taken to understand the underlying basic 
physics in a qualitative way before getting engaged in 
large scale computational efforts. Our numerical pro- 
cedure lies, hence, in the borderline of full numerical 
relativity and perturbation theory. As in p^ . the ac- 
creting matter is evolved in a curved static background 
by solving the nonlinear hydrodynamics equations. The 
response of the compact object to the infalling matter, 
which triggers the emission of gravitational radiation, is 
computed using perturbation theory. More precisely, we 
use the gauge invariant formalism of 26] and study the 
excitation of QNMs of both Schwarzschild black holes 
and neutron stars by numerically solving in the time do- 
main the even-parity perturbation equations with matter 
sources. For the case of black holes, these equations re- 
duce to the inhomogeneous Zerilli-Moncrief _3|i,.3jy equa- 
tion. One key assumption of our approach is that the 
mass of the accreting fluid is much smaller than the mass 
of the central compact object. Fluid self-gravity and ra- 
diation reaction effects are also neglected; i.e., we ignore 
the first-order metric corrections to the fluid equations of 
motion. The first approximation (no self-gravity) is in 
general valid for fluid motions in the vicinity of the com- 
pact object, where tidal forces dominate over the fluid 
self-gravity. The second approximation (no radiation re- 
action) is valid as long as the energy in the form of grav- 
itational radiation is much smaller than the kinetic or 
internal energy of the fluid. Our procedure follows then 
the same hybrid approach previously adopted in [lol ] , but 
departs from it in the formalism. We anticipate that, 
for a given numerical resolution and numerical scheme, 
the use of the Zerilli-Moncrief equation results in im- 
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proved long-term numerical stability as compared to the 
Bardeen-Press equation employed by allowing for 
an accurate computation of late time features of the GW 
signal (namely tails). The results reported in the present 
investigation are further restricted to the case of radially 
accreting shells of either dust or obeying a perfect fluid 
equation of state (EOS), where the mass density profile 
is shaped in the form of quadrupolar shells of Gaussian 
radial extent. In this respect this work can be considered 
as a necessary assessment of our numerical approach, in 
anticipation of the study of more interesting scenarios, 
namely the excitation of QNMs from perfect fluid thick 
accretion tori orbiting around compact objects, which 
will be presented elsewhere ^32] . 

The paper is organized as follows: Section Hll describes 
in some detail the theoretical framework adopted, namely 
the construction of the unperturbed stellar models, the 
general relativistic hydrodynamics equations, the pertur- 
bation equations for neutron stars and black holes, and 
the generation of time-symmetric initial data to describe 
the accreting shells. The numerical methods used for 
both the hydrodynamics equations and the perturbation 
equations are outlined in Sec. IIIII Section Hvl is devoted 
to presenting the main results of our investigation, split- 
ting the description of the black hole and neutron star 
cases into two separate subsections. Finally, Sec. sum- 
marizes the main conclusions of this work and outlines 
future directions in this research. Appendix 1X1 contains 
technical information regarding the general form of the 
source term for the inhomogeneous Zerilli-Moncrief equa- 
tion, and Appendix ^presents a comparison with point- 
like particles falling onto black holes. We use units such 
that c = G = 1. 



II. THEORETICAL FRAMEWORK 

A. Unperturbed stellar models 

The background metric of a non-rotating, spherically 
symmetric star is given by the line element 



"'^''dt^ + e^^dr^ + (rfi?^ + sin^ Mlp^) , (1) 



where a and h are functions of the radial coordinate r 
only. Assuming the star is a perfect fluid whose energy 
momentum tensor reads 



TABLE I: Frequencies of the first fluid and gravitational 
modes for the two neutron star models considered. 



Mode 


A [kHz] 


B [kHz] 


/ 


2.584 


1.666 


Pi 


3.948 


4.045 


Wl 


11.609 


10.380 


W2 


20.197 


18.429 



of hydrostatic equilibrium: 

dm 



dr 
da 
dr 
dp 



(r^ — 2mr) 



, , da 
dr=-^'+P^Tr' 



(3) 
(4) 
(5) 



where m(r) is the gravitational mass enclosed in a sphere 
of radius r. The surface of the star r = i? is determined 
by the condition p — 0. At the exterior the geometry re- 
duces to the Schwarzschild solution. The above system of 
ordinary differential equations (DDEs) can be integrated 
once the EOS is chosen and a value of the central energy 
density Cc is specified. We further assume the fluid to 
be isentropic, so that we adopt a one-parameter EOS, 
p — p{e), in polytropic form p = Ke^, where F is the 
adiabatic exponent. 

For the calculations we discuss in this work, we will 
consider two stellar models with the same gravitational 
mass, M = 1.4M© « 2.067 km, and the same adiabatic 
exponent F = 2, one model being more compact than the 
other. These two models are meant to bracket the in- 
terval of possible radii of realistic neutron stars. The 
more compact model (A) has Ec = 2.455 x 10""^^ g/cm^ 
and K = 122.25 km^, so that the radius is i? = 9.80 km. 
The less compact model (B) has ec = 0.92 x lO^^g/cm'^ 
and K = 180 km^ which leads to i? = 13.44 km. 

Our analysis is further restricted to polar perturba- 
tions induced by external matter flows. For the sake of 
reference, the first frequencies of the polar part of the 
spectrum of the two stellar models considered are listed 
in Table m They have been obtained using a frequency- 
domain code based on the classical Lindblom-Detweiler 
formulation of the perturbed Einstein equations (33] and 
described in 1341. 



B. General relativistic hydrodynamics equations 



T^j^u = [f-+p)u^Uu+pg^u, (2) 

with p denoting the pressure, e the total energy density, 
and the fluid 4-velocity, the Einstein's equations be- 
come the Tolman-Oppenheimer-Volkoff (TOV) equations 



The motion of a fluid in a curved spacetime is governed 
by the local conservation laws of baryonic number and 
energy momentum: 

J'^ = 0, V^t^'" = , (6) 



4 



where = pu^ is the mass density current and 
t^"' — phu'^u'^ + pg^'^ is the stress energy tensor for a per- 
fect fluid. In these expressions, p is the rest-mass density, 
h is the specific enthalpy, defined as h = 1 + e + p/p, and 
e is the specific internal energy. The system of equations 
is closed with an EOS p ~ p{p^ e). 

The equations of general relativistic hydrodynamics 
are a system of hyperbolic equations; as shown by |35j . 
it can be explicitly written as a system of conservation 
laws. This is accomplished by defining quantities which 
are directly measured by Eulerian observers, i.e. the rest 
mass density D — pW , the momentum density in the 
j direction Sj = phW^vj and the total energy density 
E = phW^ — p. In these definitions, W stands for the 
Lorentz factor, which satisfies W = {1 — w^)"^/^, with 
v"^ = %jv'^v^ . Here, is the 3-velocity of the fluid, de- 
fined as = /W + P"^ /a, where a and /3* are the space- 
time lapse function and shift vector, respectively, and 7ij 
are the spatial components of the spacetime metric where 
the fluid evolves. For a generic spacetime, the system of 
equations we solve reads [s^ 



\ dt dx^ 



= S(w), (7) 



where g — det(g^i,), U(w) = {D,Sj,E — D), and 
w = {p,Vi,e) is the vector of primitive variables. The 
expressions for the flux and source vectors, F'(w) and 
S(w), can be found in Ref. [s^. We note that in the cur- 
rent work we use the standard form of the Schwarzschild 
metric in polar-radial coordinates (t, r, t?, ip) to describe 
the exterior spacetime of the compact object. Hence, the 
above expressions are specialized accordingly. 



C. Stellar polar metric perturbations induced by 
hydrodynamical sources 

1. Interior equations 

A formulation of the equations describing the polar 
perturbations of a star in the Regge- Wheeler gauge [s^ . 
written in a form suitable for numerical simulations in 
the time domain, was first discussed by . An alterna- 
tive derivation of the same equations, based on the lin- 
earized Arnowitt-Deser-Misner (ADM) formalism, can be 



found in Recently, building on the work of Gerlach 
and Sengupta Martfn-Garcfa and Gundlach have 

developed a gauge-invariant and coordinate-independent 
formalism for non-spherical perturbations of spherically 
symmetric spacetimes .25..^ 2 6 ,.27J. In our work, we follow 
the formalism laid out in [2y| specifying the equations to 
the case of a static spherical star and choosing the Regge- 
Wheeler gauge The equations we obtain are then 

equivalent to those of Refs. and 0|, although differ- 
ent metric variables are used. Hence, for each {I, m) pair, 
the even-parity metric perturbation Sg^^i, is parametrized 
by two scalar quantities, k (the perturbed 3-conformal 
factor) and x (the actual gravitational wave degree of 
freedom), so that it reads 



fix + fc)e2' 


V 















kr^ sin^ i9 / 



(8) 



Since the background is static, "0 is not an independent 
quantity, as it can be obtained by quadrature from k and 
X H^l; the relationship with the usual Regge- Wheeler 
variables can also be found in We note that, al- 

though the polar problem on a static star is known to 
have only two metric degrees of freedom inside the star 
and one degree of freedom outside (represented by the 
Zerilli-Moncrief function), we have decided to consider 
an additional variable inside the star, the perturbation 
of the relativistic enthalpy H = dp/(p + e), as suggested 
by earlier studies of the subject [ifil ll^ . Correspond- 
ingly, at the exterior we evolve two (constrained) degrees 
of freedom instead of just one (see below). We notice, 
however, that successful evolution algorithms using the 
actual number of degrees of freedom have been developed 
in the past in more general frameworks |2^. 

We formulate, then, the polar perturbation problem 
through a couple of hyperbolic equations for x and H plus 
an elliptic equation, the Hamiltonian constraint, which is 
used to compute k at every temporal slice. This permits 
us to obtain the frequencies of the stellar pulsation modes 
with an accuracy comparable to frequency domain calcu- 
lations. The set of equations reads 
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-e^H -2 



2e 



2b 



2 ^ 6m 
+ 87re - — 



(x + fc) H 5— X 



2 ?Ti 

47rr(5p - e) h 10— 



(9) 



m , o 

;:2 (1 + ^? 



47rrp (1 - 2cl) + ^iirre - 



2 (m + Anpr^y 
(r — 2m) 



47r (3p + e) 



(x + fc) 



(10) 



-2b I 



]k + 











-26 



3 m 



87r(p- 



-H = 0, 



(11) 



where X = l{l + 1) and = dp/de is the sound speed. 
At the star surface e, p and vanish, and m(i?) = M; 
thereby the evolution equation for H reduces to the ODE 



M{R - 2M) 



2M^ 



(x + fc) .(12) 



Since we are setting F = 2, the term proportional to H 
in the Hamiltonian constraint l)ll|l is regular at r — R. 



monies reads 



00 I 

^M'' = X! X! 

00 / 



(13) 



EE 



Exterior equations with a general source term 

The exterior equations in vacuum are readily obtained 
from the interior equations setting p = p = and 
m = M. As we have seen, the equations for the in- 
duced metric perturbations are basically wave equations 
with potentials. The presence of an extended object out- 
side the star or black hole reflects in the fact that these 
equations are not homogeneous anymore, but they con- 
tain source terms involving the stress-energy tensor of 
the external fluid, t^j^. Whereas for point-like particles 
these source terms can be explicitly computed by ana- 
lytic techniques 0,1131, in the case of extended objects 
they involve steps that cannot be managed analytically 
[Toj . To the best of our knowledge, an explicit expression 
of the source terms for a general stress energy tensor has 
not yet been reported in the literature. 

Following the notation of Ref. [2^, the (gauge- 
invariant) decomposition of t^^ in polar spherical har- 



where the capital indexes run over the M"^ Lorentzian 
manifold and the lowercase indexes over the unit radius 
2-sphere S'^, as the background spacetime can be written 
as a direct product M"^ x S'^. We also follow for the 
deflnition of the scalar and vector spherical harmonics, 
y£m respectively, and of the tensor spher- 

ical harmonic, Zf'j," = K^'^" + ^jabY'^"'. Here, the nota- 
tion : a stands for the covariant derivative with respect 
to the metric ^ab = diag(l, sin^ i9) of 5^. In Ref. |25l| . 
the homogeneous equations of linearized nonspherical po- 
lar perturbations of a general time-dependent spherically 
symmetric spacetime were obtained. If we choose the 
background spacetime to be the Schwarzschild solution 
and interpret tfj_i, as a source induced by a certain dis- 
tribution of matter on this spacetime, one can arrive af- 
ter some algebra at the polar perturbation equations of 
a Schwarzschild spacetime with source terms, which are 
given by 
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r 2M 


2e-^ - 6 


1 ~ — 


r 



,x -^-2 2 /5M , 
{x + k) + -^X - 1 ) X.r 



(14) 



2 3M\ , e-^fc 

T 

r r'^ r 



I 



(15) 



The sources S-^, and S'-h read 



e2("-'') J r/r + (rj™) - 2 (rf") ^ + - ( lll^e^" - 3 I (t; 



C . ^ — 2arTiim 



1 /5M^2, 



1 



\ r 2 



o2b 



rphn ^2brj-^im 



(16) 
(17) 



It is known that the perturbations of the Schwarzschild 
spacetime are described by a hyperbohc equation for a 
single function, originally written in the frequency do- 
main and in the Regge- Wheeler gauge by Zerilli 30] and 
later in the time domain by Moncrief |3l| . but adopting 
a gauge invariant formulation. Following the normaliza- 
tion convention of 17], the Zerilli-Moncrief function is 
related to k and x as follows ^3 



4r^e 



2„-2fc 



A [(A - 2)r + 6M] 



2 r 



(18) 



In presence of matter sources, this function is a solution 
of the inhomogeneous Zerilli-Moncrief equation 



Z.tt ~ e^^'^-^^Z.,, = ^e^^Z, + ViZ + S, 



(19) 



where the Zerilli potential V/ is given by 



Vi 



(20) 



A(A - 2Yr^ + 6(A - 2YMr^ + 36(A - i)K'Pr + 72M' 
r3 [(A - 2)r -t- &Mf 



The source term Sz appearing in Eq. (|19l) for a general 
matter distribution has not been reported in the litera- 
ture in a form suitable for time domain computations. Its 
derivation can be found in Appendix A. The final result 
is 



Sz - 



A [(A - 2)r + 6M] 1 (A - 2)r + 6M 



00 



,-26 



2Mr + r2(A - 4)1 rf{" + 2r^ {T^^) ^ - 2r{r - 2Mf (t/]") ^ 



4A(r - 2M)Tf 



2A(1-^Va^ 

r 



r / 



Tl" + 4r(r - 2M)T^ 



(21) 



In the case of a point particle, this general source term reduces to that given in p^. 
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For any multipole {l,m) the radiated energy is com- 
puted from the Zcrilli function. It is given by 



E 



— dt 
dt J 



647r {1-2)1 J_^\ I ^ ^ 

By defining the Fourier transform Z of Z = Z^"^ as 

/oo 
e-''^'Z{t,r)dt , (23) 
-oo 

we obtain the energy as follows: 



dio 



1 (/ + 2)! 



647r2 (l~2)\ Jo 



duj 



(24) 



D. 



Black hole polar metric perturbations induced 
by hydrodynamical sources 



When the neutron star is replaced by a Schwarzschild 
black hole, the polar perturbation problem becomes 
much simpler, since one only needs to solve the inho- 
mogeneous Zerilli-Moncrief equation. As mentioned in 
the Introduction, in order to draw a comparison with the 
neutron star ca.se, we consider the black hole case dis- 
cussed in Ref. [ij, but solving Eq. (|19|l instead of the 
inhomogeneous Bardeen-Press (BP) equation. The un- 
derlying motivation behind this choice is the possibility 
of performing long-term stable evolutions that allow for 
the extraction of late time features (radiative power-law 
tails) in the GW signals. 

Working with the same numerical method, this result 
seems to be unreachable with the BP equation, because 
this equation is intrinsically unstable. To make the ar- 
gument clearer, let us recall that this equation, written 
using the Regge- Wheeler tortoise coordinate r* [s^ . 



r* = r + 2Mlog (— 1 



). 



(25) 



reads 



4(r - 3M) 



A 



6M 



+ (1 + 2) {l-l) 



Y + — r , 

r 

(26) 



where the (complex) function Y is related to the Weyl 
tensor tetrad component by Y = A is the horizon 
function A = — 2Mr and T is the source term deter- 
mined by matter flows jlOj| . The polar metric perturba- 
tions correspond to the real part of Y, the axial ones to 



its imaginary part. At r = 3M the term proportional 
to Yt changes sign: depending on the relative sign be- 
tween Y^tt and Ft , one may view Ft either as a damping 
term (when the signs of both coefficients agree, i.e. for 
r < 3M) or an antidamping term (otherwise). As argued 
by Krivan et al. 37], who analyzed in detail the general 
case of the Teukolsky equation for Kerr black holes, this 
is likely the origin of exponentially growing modes that 
appear when the equation is numerically solved by stan- 
dard finite-differencing explicit methods. Although some 
attempts to delay in time the onset of the instability have 
been investigated in the literature [lol l37j , it remains an 
open issue. As noted above, this effect is particularly 
disturbing in the presence of matter sources, since the 
instability is always occurring before that the late-time 
state of the system is reached. 

On the other hand, using the tortoise coordinate, 
Eq. lO reads 



ViZ + S, 



(27) 



where the source term is the same one given by Eq. (|21ll , 
with the derivatives with respect to r consistently re- 
placed by derivatives with respect to r*. This equa- 
tion does not present terms which may cause exponential 
growing modes, and, therefore, it permits stable evolu- 
tions. 



E. Initial data 

The nontrivial issue of how to specify suitable initial 
data (i.e. gravitational radiation free) in the presence of 
sources has been addressed to some extent in a number 
of works [H [ll HI m Sf. In particular, initial data 
suitable to describe point-like particles scattered by stars 
or falling onto black holes can be found in and [s^ . 
respectively. The common procedure is to choose ini- 
tial data such that the Hamiltonian and momentum con- 
straints are satisfied at the initial temporal hypcrsurface. 
If the matter source is initially at rest, the initial condi- 
tions are time symmetric and the momentum constraint 
is automatically satisfied if the Hamiltonian constraint 
is. However, if velocity fields are present initially, the 
momentum constraint must be solved for too [T3| . 

In the simulations reported in this work, we choose 
time-symmetric initial configurations where the matter 
distribution, shaped in the form of quadrupolar shells of 
dust or perfect fluid, is falling onto the central compact 
object from rest. Hence, we only need to consider the 
Hamiltonian constraint, Eqs. Hll|l and IjlSfl . This con- 
straint is a single equation for three unknowns, H, fc, 
and X- -By setting H = 0, one of the functions k or x 
can be specified freely and the constraint is then solved 
for the remaining one. Furthermore, the initial data are 
chosen so as to minimize the amount of gravitational ra- 
diation present initially. This is done by choosing x — ^ 
and solving for k with a given source Tqq |38j . Given fc, 
the initial profile of Z is then computed using Eq. H18|l . 
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While this prescription should ensure that the initial data 
are free of any spurious GW content other than that due 
to the presence of the matter source, in practice this is 
not exactly the case, a certain amount of GWs being al- 
ways present. Its origin is related to the finite value of 
the initial location of the shell. For such a configuration, 
the impossibility of specifying the GW contribution asso- 
ciated with the shell in a way consistent with its past his- 
tory provokes a transient burst in which the excess GWs 
are radiated away. Similar situations were considered in 
Refs. Elllllil, 

where the source of the perturbations 
was a particle orbiting around or scattered off a star or 
a black hole. In order to minimize this problem in our 
simulations, we freely evolve the perturbations without 
the hydrodynamics part until the initial unphysical burst 
of gravitational radiation leaves the numerical domain, 
taking the final profiles as the initial state for the actual 
simulations. In the black hole case, this simple procedure 
permits us to avoid completely any kind of initial data 
interference. For stars, however, this approach triggers 
the oscillations of the fluid modes. Hence, we further pro- 
ceed by resetting 7J = 0, in order to obtain the correct 
initial model. 



III. NUMERICAL FRAMEWORK 

The numerical schemes we have implemented to solve 
the hydrodynamics and perturbation equations (outlined 
in the preceding section) are used with some technical dif- 
ferences in both scenarios under study, neutron stars and 
black holes. The first difference is that, in the neutron 
star case, the overall grid is uniformly spaced in the ra- 
dial coordinate r, while for the black hole it is uniformly 
spaced in the tortoise coordinate r* . The numerical do- 
main chosen to discretize the hydrodynamics equations is 
always smaller than that of the perturbation equations, 
which is extended in both directions, towards the hori- 
zon of the black hole (or the origin of coordinates in the 
case of the star) and towards large radii. This procedure 
avoids or minimizes the effect of the spurious reflection of 
waves at both boundaries. For stars, the hydrodynamics 
domain begins at the first cell outside the stellar surface, 



and extends up to 



= 108 km. We choose the same 



resolution (Ar ^ 0.07 km) for both stellar models, so 
that the interior is covered with 146 points for model 
A and 200 points for model B which, we recall, is less 
compact. The hydrodynamics grid is covered with 1400 
zones, a resolution chosen to ensure convergence. The 
external wave domain extends up to r ~ 1500 km and is 
covered by roughly 22000 zones. On the other hand, in 
the case of black holes the hydrodynamics domain starts 
very close to the horizon so as to include as much as 
possible the peak of the potential barrier as well as its 
falloff toward r* = — oo. The sensitivity of our numeri- 
cal results to the location of the inner boundary of the 
hydrodynamics domain is discussed in Sec. II VI below. 
Before turning to describe the numerical schemes im- 



plemented in the code, it is worth commenting that in a 
realistic scenario the mass of the compact object would 
grow in time as the accretion process proceeds. Such 
a possibility, however, has not been considered in the 
present simulations. Technically, this would require re- 
computing at each evolution time step the equilibrium 
structure of the star with a mass M + 5M{t) (or, anal- 
ogously, increasing the black hole mass). In our simula- 
tions we simply assume that 6M ^ M, neglecting the 
effect of a dynamically growing mass. 



A. Evolution of the external fluid 

The hydrodynamics part of our code is the same as 
that used in the simulations reported in Ref. [lOf . In this 
code in order to evolve dynamically the infalling fluid 
shells, the general relativistic hydrodynamics equations 
are solved using a Godunov-type scheme (see e.g. [s^ 
for details). In axisymmetry [d^ = 0), the vector U of 
evolved quantities appearing in Eq. ((7|) is updated from 
time level to i""*"^ according to a conservative algo- 
rithm 



U"+i = U" . - — 



At 



-l/2j 



-l/2,i 

f AiS 



(28) 



where At = — t", and Ar and AS indicate the radial 
and angular grid spacing, respectively. In practice, a con- 
servative, second-order, two-step Runge-Kutta algorithm 
is employed instead of Eq. H28(l . In the above equation 
i and j label the radial and angular zones, respectively. 
The numerical fluxes (e.g. F^_^.i/2 j) are calculated at ev- 
ery cell interface using an approximate (linearized) Rie- 
mann solver built upon the characteristic information of 
the Jacobian matrices of the system. The reader is ad- 
dressed to '3^ for further details. 

The matter model we choose for the accreting shells 
can be either dust or perfect fluid. In the former case, 
as p = 0, the Riemann solver must be specialized ac- 
cordingly to avoid divergences. For shells accreting onto 
black holes we employ dust, not only to allow for a di- 
rect comparison with the results of [13 , but also because 
of numerical difficulties encountered when using a per- 
fect fluid in regions very close to the event horizon. In 
Schwarzschild coordinates some metric components blow 
up at the horizon, which affects the evolution of hydro- 
dynamical quantities. In particular, the coordinate flow 
velocity becomes ultra-relativistic and reaches the speed 
of light at the horizon, making the Lorentz factor infinite. 
If the inner boundary needs to be placed very close to the 
horizon to capture the fall off of the Zerilli potential (we 
have checked that r* = —BOM is a reasonable value; see 
below) the metric and hydrodynamical quantities, de- 
spite being regular, present steep radial gradients, which 
make the numerical evolution difficult. In our current 
code and with the maximum grid resolutions affordable. 
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we have found severe limitations in the perfect fluid ac- 
cretion case to move the innermost boundary past values 
of about r* — 3M which, as we show below for the case 
of dust shells, are not yet close enough to the horizon to 
capture the gravitational wave emission unambiguously. 
We note that this numerical drawback can be entirely re- 
moved by using horizon- adapted coordinate systems such 
as ingoing Eddington-Finkelstein coordinates employed 
in simulations of perfect fluid accretion onto black holes 
in Ref. Such a procedure for perfect fluid accretion 
may be attempted in future work. 

However, for the case of neutron stars we do not en- 
counter numerical difficulties in evolving perfect fluid 
shells. Therefore, in those simulations we do not consider 
quadrupolar dust shells, but only perfect fluid shells. 
This is also further motivated from the reflecting bound- 
ary conditions we impose at the surface of the star (see 
next paragraph), which may naturally lead to the ap- 
pearance of shock waves which could not be treated in 
the case of dust. 

The boundary conditions to impose upon the hydro- 
dynamics variables completely depend on the case under 
study. In the black hole case, we adopt ingoing radial 
boundary conditions at the innermost grid point of the 
hydrodynamics domain, chosen as close as possible to 
the event horizon. As for the neutron star, the boundary 
conditions to impose to the infalling fluid shell as it ar- 
rives at the stellar surface should take into account the 
interaction between the shell and the star. As a result 
of the complexity in modeling this phenomenon within 
our perturbative approach, we impose reflecting bound- 
ary conditions at the inner edge of the hydrodynamics 
domain (i.e. the surface of the star), so that the star is 
seen by the external fluid as a hard sphere. This choice 
includes the most relevant effect; that is, the pressure 
gradient stops the infalling matter and avoids the viola- 
tion of energy conservation that outgoing boundary con- 
ditions would introduce (as happens in the black hole 
case when the inner boundary is not close enough to the 
event horizon) . Correspondingly, at the outermost radial 
boundary we impose stationary values of the so-called 
Michel solution at those grid points ^3 ^-t times. 
We note that the stationary Michel solution is used in 
the entire hydrodynamics grid in order to provide a dy- 
namically unimportant, spherically symmetric accreting 
atmosphere surrounding the fluid shell. In the angular di- 
rection, axial symmetry fixes the appropriate boundary 
conditions at the axis (i) ~ and tt). 



B. Integration of the perturbation equations 

1. Black hole case 

We can solve the Zerilli-Moncrief equation using either 
a standard three-level leapfrog method, at second-order 
in both space and time, or a second-order Lax-Wendroff 
scheme. Both methods have proved robust and stable 



enough for long time evolutions. In particular, they al- 
low the computation of long-term features in the GW sig- 
nal, namely the distinctive power-law tails following the 
black hole ringdown. However, when the leapfrog method 
is used, we find some high-frequency numerical noise of 
small amplitude at the very end of the tail. This numer- 
ical noise is not present when using the Lax-Wendroff 
scheme. Therefore, all results reported below for black 
hole spacetimes are obtained using this latter method. 
We note, however, that the second-order leapfrog pro- 
duces noise-free results in the neutron star case. The 
particular form of such a scheme is discussed in the fol- 
lowing section. 

In order to apply the Lax-Wendroff method, the 
Zerilli-Moncrief equation is written as a first order hy- 
perbolic system as follows; 



dtV + dr'F = S, 



(29) 



with 



U 





where we have introduced the variable w = Z^t 
The time update algorithm is given by 



^ 2Ar* 



At' 



2Ar2 



-j+i 



2F" 



where the matrix A reads 




AiS" , (31) 



(32) 



To ensure stability of the code, the time step At must sat- 
isfy the usual Courant-Fredrichs-Levy (CFL) condition, 
i.e. At = AcFbAr*, with Acfl < 1- 

The same radial resolution is used for the two domains 
(hydrodynamics and perturbations) present in the com- 
putational grid, namely Ar* ~ 0.04M. The angular do- 
main of the hydrodynamics grid extends from to tt 
and it is covered with 20 zones. As mentioned before, 
the hydrodynamics domain is within the wave domain, 
which is much larger. It extends from r* ~ — 50M to 
r* — 30M and is covered with about 2000 zones. The 
peak of the Zerilli potential is located at r ~ 3.1M 
(r* ~ 1.9M). Since the potential decays exponentially 
towards the horizon, r* = —BOM is a sufficiently small 
value to minimize the effects of the truncation on the 
GW signal (see below). On the other hand, in the pos- 
itive r* direction it is not necessary to extend the hy- 
drodynamics grid much beyond the position of the cen- 
ter of the shell, because the shell is always collapsing 
towards the black hole. Correspondingly, the grid for 
the perturbation equations extends from r* ^ — 876M 
to r* ^ 1250M. Standard Sommerfeld outgoing-wave 
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conditions are imposed at the external boundary, and 
ingoing-wave conditions at the black hole horizon. Ob- 
taining an optimal resolution in the simulations requires 
us to use 2 X 10^ points from r* — towards the horizon 
and some 3 x lO'* zones in the opposite direction. Such 
a large number of zones, however, does not imply any 
numerical limitation in the code as solving the Zerilli- 
Moncrief equation is a one-dimensional problem. 

2. Neutron star case 

The numerical algorithm used for solving the stellar 
perturbation equations does not use the gauge-invariant 
metric perturbation variable X: but instead S = x/^ 
(the equations are rewritten accordingly). The reason for 
this modification is that x grows proportionally to r for 
r — > oo, while the amplitude of S remains finite. This is 
an important property that helps to avoid unphysical in- 
stabilities of a numerical nature. The evolution equations 
for the independent variables S and H are discretized on 
a uniformly spaced grid and solved using the three-level 
leapfrog method. The remaining elliptic equation for k is 
also discretized in space and then written as a tridiagonal 
linear system, which is solved at each time step inverting 
the corresponding matrix by a standard LU decomposi- 
tion. 

In order to avoid numerical problems at the origin, the 
left interface of the first grid zone is chosen to coincide 
with the center of the star, r — 0, while the surface of 
the star, r ~ R, is located at a cell center. The first 
point of the grid, j = 1, is then located at r = Ar/2, 
where Ar is the radial grid spacing. Correspondingly, 
the surface, which is labeled by the J cell, is located at 
R = Ar/2 + {J—l)Ar. As mentioned before, S is evolved 
both inside and outside the star, while H is only evolved 
inside the star according to Eq. (|10|l and at the surface 
with Eq. H12(l . The same grid spacing is chosen inside 
and outside the star. 

The radial derivatives of S, k, and H are discretized 
to second order, 

[A.rr]; = (A^+i - 2A- + /Ar\ (33) 

[Ar]] - {A]+, - A-_,) /2Ar (34) 

(with A = S,k OT H), but at j — J we use backward 
second-order differencing for H,r as 

[HX = {Hj^2 - ^Hj^i + iHj) /2Ar . (35) 

Similarly, the second time derivatives of S (and H) are 
approximated as 

[S.ttT^ = (5;+' - 25; + S]'-^) /At\ (36) 

As for the black hole case, we impose standard Som- 
merfeld outgoing-wave boundary conditions at the out- 
ermost radial zone. Partial reflection from the external 
boundary is still present, but it can be minimized (or 



causally disconnect its influence) by placing the outer- 
most point sufficiently far. At the origin of the radial 
coordinate r = 0, all fields are regular and vanish [26l |. 
For Z = 2, as the origin is chosen to lie at the interface of 
the first cell, there is no need for regularizing the fields, 
as had to be done in 0, 0| . 

In order to have an internal consistency test, we have 
studied the evolution of a generic enthalpy profile in the 
source-free case, computing the Zerilli-Moncrief function 
in two different ways. In the first option, it is computed 
algebraically at every temporal slice from k and x by 
Eq. H18|) , all over the external domain; in the second pos- 
sibility, it is matched to k and x using Eq. (|18|l just at 
the stellar surface, and then it is evolved independently 
in the exterior region by solving Eq. (|19|) . Both solutions 
agree very well, and only minor differences are observed 
during transient states. 



IV. RESULTS 
A. Black hole simulations 

As mentioned in the Introduction, time-dependent nu- 
merical simulations of the accretion of dust shells falling 
isotropically onto a Schwarzschild black hole were pre- 
sented in [lOl. The aim of those simulations was to char- 
acterize and estimate the gravitational radiation emitted 
in accretion events of finite-size objects and to compare 
the outcome with the point particle case amenable to 
semi-analytic investigations. 

The time domain simulations discussed in were 
based on the Newman-Penrose formalism, from which a 
master wave equation for black hole perturbations was 
derived by Teukolsky [ 1 1| . As mentioned above, in the 
case of a Schwarzschild black hole, the Teukolsky equa- 
tion reduces to the BP equation. Here, we start by 
reexamining the GW emission of a Schwarzschild black 
hole as a result of the radial accretion of a quadrupolar 
dust shell. In this work, however, we use the inhomoge- 
neous Zerilli-Moncrief equation for two different reasons: 
firstly, because it provides a direct comparison with the 
neutron star case and, at the same time, allows for cross- 
checking our results with those of ; secondly, the main 
motivation in reexamining this problem was the fact that 
the Zerilli-Moncrief equation allows for long and stable 
time evolutions. This property permits the investigation 
of long time features in the GW spectrum. 

We consider an external matter distribution of total 
mass /i, much smaller than that of the black hole, namely 
/i = O.OlAf. As in 10], the shell density distribution is 
parametrized according to 

p = Po+Pm..e-'^^--'-^'^\mH, (37) 

where rp is the initial position of the center of the shell 
and K, controls its width. The background density po is 
chosen to be very small (~ 10^^^ km ~^) to simulate the 



11 




20 40 60 80 100 o 50 100 150 200 

u/(2M) u/(2M) 



FIG. 1: Gravitational waves emitted by a Schwarzschild black hole excited by infalling quadrupolar shells of given compactness 
K [see Eq. (I37ll 1 . The right panels show the logarithm of the GW signals shown in the left panel in order to highlight the onset 
of the fundamental QNM ringing. The shells are falling from ro — 15M and the observer is located at robs ~ 125Af from the 
origin of coordinates. 



vacuum outside the black hole. The angular structure 
of the matter density is assumed to have a quadrupo- 
lar profile, being given by sin^-!^. The value of pmax is 
obtained from the condition that the volume integral of 
Eq. (|37() equals /i. Papadopoulos and Font analyzed 
the dependence of the black hole QNMs excitation on 
the various parameters of the shell, namely its width, 
initial location, and its initial velocity field. Here, we be- 
gin by studying the excitation of the black hole QNMs 
varying the compactness of the shell, i.e. its width k 
and the position where it is at rest. We restrict our- 
selves first to shells which are falling from a fixed location 
ro = 15M and are initially at rest. The distant observer 
is located at robs — 125 A/ from the origin of coordinates 
(r* - 133. 5A/). 

The response of the black hole as a result of the ac- 
cretion process when shells of different widths k are con- 
sidered is shown in the left panel of Fig. ^ in which we 
plot the time evolution of the Zerilli function. Note that 
in this figure we use retarded (observer) time u = t ^ r* . 
The right panel shows the absolute values of the wave- 
forms, but in logarithmic scale and separating the dif- 
ferent models for clarity. The values of k are chosen in 
order to cover the range in which qualitative differences 
among the signals are observed. The process can be di- 
vided into three phases: (i) motion of the shell before 
the bulk reaches the peak of the Zerilli potential, from 
u/2M ^ —10 up to ujIM ~ 50, where the emission is 
purely due to the variation of the quadrupolar moment 
of the shell; (ii) shortly after u/2M ~ 50, when the bulk 
of the shell approaches the peak of the potential, Z ex- 
periments a rapid variation due to the interaction with 
the potential barrier, resulting in a burst-like short sig- 



nal. For more compact shells, higher amplitudes of the 
Zerilli function are obtained, (iii) The last part of the 
signal, from u/2M > 55 — 60 onward, is characterized 
by highly damped oscillations at early times (the funda- 
mental QNM ringing of the black hole), followed by the 
power-law radiative tail at the end. In the right panel 
of Fig. n it is visible how the onset of the fundamental 
QNM ringing occurs only when k > 0.3; i.e., wider (less 
compact) shells do not succeed in exciting the fundamen- 
tal mode of the central black hole, but rather the signals 
show a large wavelength oscillation. This feature is in ex- 
cellent agreement with the results reported by 'T0| using 
the BP equation. If k is too small, the impinging GWs 
cannot be fully transmitted beyond the potential barrier 
to cause the ringing of the fundamental QNM. Thus, the 
large wavelength oscillation is determined by the gravi- 
tational pulse driven by the shell, which is almost com- 
pletely reflected back by the potential barrier. It is worth 
mentioning that, also in this case, the black hole space- 
time reacts to the external perturbation, as we found that 
the other QNMs, characterized by frequencies lower than 
the fundamental one can be slightly excited during 
the process. However, as a result of their high damping 
times, the energy that can be released through GWs is 
negligible with respect to that generated by the variation 
of the quadrupolar moment of the shell during the infall 
(roughly three orders of magnitude lower for k — 0.1). 
We note that the interaction of the shell with the black 
hole potential, even for small values of k, is confirmed by 
the presence of the distinctive late-time decay. 

It is known that the generation of the fundamental 
QNM ringing is associated with the peak of the poten- 
tial The bulk of the shell crosses the peak at time 
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FIG. 2: Energy distribution for some selected values of k of 
the waveforms shown in Fig. Q The power spectrum has a 
complex shape characterized by various distinctive peaks re- 
sulting from the interference between emitted and backscat- 
tered GWs (see text for details). 



t/2M ~ 50, which is at retarded time u/2M ~ 49. As 
argued before, the maximum GW amplitude and the fol- 
lowing QNM ringing are reached somewhere around this 
point. The end of the accretion process occurs at time 
tjlM ~ 79.5, when the center of the bulk tq reaches 
the innermost boundary of the hydrodynamical numer- 
ical domain. This time corresponds to a retarded time 
m/2M- 104.5. 

The energy spectra for some selected values of k are 
shown in Fig. [21 For all values of k, the spectrum dis- 
plays a complex structure, with several distinctive peaks 
that are more or less evenly spaced. The spacing between 
the maxima is roughly given by 0.1 (in units of IMlS). 
The presence of these peaks is interpreted as an effect of 
the interference between the gravitational waves emitted 
by the shell during its motion and the radiation emitted 
at earlier times, which has already been backscattered 
by the potential. This interference fringes are not new. 
In fact, similar patterns were found by Lousto and Price 
|43| when studying the signal emitted by a point parti- 
cle falling radially onto a Schwarzschild black hole from 
a finite distance tq. However, some differences between 
their case and ours must be stressed. In Ref. "Zs^ , the ini- 
tial data had some GWs content, and the authors argued 
that the evenly spaced bumps found in the energy spec- 
tra (whose amplitude and spacing depended on rg) were 
mainly due to the interference between the initial data 
pulse and the GWs emitted by the infalling point parti- 
cle during its motion. In Appendix ^ we confirm this 
result using a time domain code, by comparing the GWs 
energy spectra generated by a falling point particle with 
and without the initial GW content. In the latter case, 



FIG. 3: Excitation of the I — 2 fundamental QNM as a func- 
tion of K. The spectrum is computed selecting only the part 
of the signal corresponding to the ringing. The dashed verti- 
cal line indicates the fundamental QNM frequency {2Mlu « 
0.7473) of a Schwarzschild black hole 0. Note that for 
K — 0.1 the spectrum is barely visible at the bottom of the 
plot. 



we find that the amplitude of the interference bumps is 
reduced. For the case of an imploding shell, despite the 
initial data effect being eliminated by construction, we 
find well-defined interference fringes in the energy spec- 
tra. The extended size of the matter distribution results 
in the existence of interference patterns even when the 
GW initial content is minimized. 

As for the case of point particles ^3 j notice that the 
part of the waveforms which strongly contributes to these 
features in the spectrum is that extending up to the sec- 
ond zero of the signal {u/2M ~ 42), just before the burst. 
This portion of the waveform carries the imprint of the 
radiation emitted during the accretion process, before the 
bulk of the shell crosses the peak of the Zerilli potential. 
This is confirmed by the results shown in the spectra of 
Fig. 121 obtained considering only that part of the signal 
from the third zero of Z onwards. This helps eliminate 
as much as possible the GWs contribution related to the 
motion of the shell and to dig out the actual QNM ring- 
ing signal. The spectrum of this late time signal closely 
corresponds to that of a Schwarzschild black hole radi- 
ating via the fundamental mode {2Mlo k, 0.7473), whose 
frequency is marked by the vertical dashed line in the fig- 
ure. The width of the peak is consistent with the damp- 
ing time of the fundamental mode (see ^3), confirming 
that, if the black hole is successfully excited, the energy 
is radiated mostly in the fundamental mode. 

In Fig. 21 we show the superposition of the QNM ring- 
downs for some selected values of k. As expected, the 
power law of the late time tail does not depend on the 
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FIG. 4: Black hole QNMs ringdown and tails for various shell widths k. The left panel shows that the GW burst resulting 
from the bulk of the shell crossing the peak of the Zerilli potential occurs at retarded time u/2M ^ 49.5 (first vertical dashed 
line). The late-time power-law tails are all perfectly superposed for the various values of k considered. At u/2M ^ 104.5 
(second vertical dashed line), when the shell leaves the hydrodynamical numerical domain, a second ringing appears. This 
ringing is a numerical artifact (see text for details). The right panel depicts on a log-log plot the late time behavior of the 
longest simulation for a shell with k = 10. A fit to the solid line gives Z ~ (tj'lM^'^'' ^ in excellent agreement with the analytic 
late-time faUoff derived by Price [i^, Eq. 



shell width /t, being all nicely overlapped. This important 
feature was not accessible to the simulations reported 
in [To| , due to the appearance of numerical instabilities 
when solving the BP equation for sufficiently long evolu- 
tion times. The late time behavior of gravitational per- 
turbations was first studied in detail by Price '45'] using 
analytic techniques. The most recent and exhaustive dis- 
cussion of this topic can be found in . The gravita- 
tional multipole perturbations with I > 2 are expected 
to fall off at large t as 

/ ^ -{21+3} 

i.e. as ~ in our case {I — 2,m — 0). The best fit to the 
tails shown in the left panel of Fig. 01 which correspond 
to simulations that extend up to ~ 3.67 ms of evolution 
{u/2M ~ 209), gives Z ~ {t/2M)-'^-^^. The somewhat 
large difference with the expected analytic value given 
by Eq. H38|l is due to the fact that the signal has not 
reached yet the late-time state. In order to prove this, 
we perform a much longer run (up to u/2M 609) em- 
ploying a larger grid of 7 x 10* zones, which corresponds 
to roughly twice the number of zones used in the previ- 
ous simulations. For this run, we set n = 10, keeping 
the same values for the remaining parameters. The right 
panel of Fig. 0] shows the results of this long run. The 
best fit to the late time signal is now Z ~ (t/2M)~^ °^, 
in close agreement with the expected value. In order to 
further improve the numerical results it simply suffices to 
perform even longer simulations employing larger grids. 



Some more comments are relevant about Fig. ^ For 
K = 1, 10, and 100 the signals present a strange feature, 
with a second ringing starting at u/2M ~ 104.5 (indi- 
cated by a thick vertical dashed line in the left plot). We 
have checked that this second ringing starts at the time 
when the center of the shell leaves the numerical domain 
through the innermost boundary. The appearance of this 
ringing seems to be a numerical artifact, as a small dis- 
continuity in the fields is introduced when the center of 
the shell leaves the grid. The black hole reacts to this by 
emitting GWs in the form of the second ringing, until the 
late time tail is reached. This unphysical ringing could be 
avoided by extending the numerical domain as much as 
possible towards r* — —oo (i.e. towards the event hori- 
zon). In practice, since the accretion processes we are 
interested in happen outside the horizon, where the peak 
of the Zerilli potential stands, it suffices to choose the 
innermost boundary of the hydrodynamics domain such 
that the exponential falloff of the potential is properly 
captured. The observation that most of the energy is 
released at low frequencies shows that the choice of the 
radial extent of the hydrodynamics domain with respect 
to the width of the Zerilli potential is of paramount im- 
portance to obtain the correct GW signals and the corre- 
sponding energy spectra. This effect is more important 
in the case of an extended shell, whose size changes with 
time due to the presence of tidal forces which tend to 
disrupt it before being swallowed by the black hole. It 
is the complex interaction with all the structure of the 
potential which determines the GW emission. 
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FIG. 5: Dependence of the GW energy spectra on the radial 
location of the innermost boundary r'^i^ of the hydrodynamics 
numerical domain. The spectra plotted correspond to k = 10 
quadrupolar shells falling from ro = 15M, the GWs being 
extracted at robs ~ 125Af. The unphysical high-frequency 
components disappear as r^j^ becomes smaller than about 
— 20M. The spectra converge for values of r^in as small as 
-50M (solid line). 



In order to study how the radial extent of the hydrody- 
namical domain affects the waveforms and energy emis- 
sion, we focus on an accreting shell of fixed width k = 10 
and vary the radial location of the innermost bound- 
ary The results of these simulations are shown 
in Fig. We start with r'^-^ = —3M and gradually 
push r*jjjj towards the event horizon, choosing the val- 
ues -lOAf, -20M, -40M, and -50M. Although the 
low-frequency part of the spectra {2Mlu < 0.4) is al- 
most unaffected by the location of the inner boundary, 
high-frequency components become evident when r^^^^ is 
larger than — 20M. The energy spectra for r*^i^ = — 40Af 
and r*j^;„ = —50M are in practice identical. This means 
that such a radial extent suffices to capture all the rel- 
evant GW physics of the accretion process. However, 
using a less conservative value for the radial location of 
the inner boundary, e.g. r^^^ = — 3M, results in an arti- 
ficial spectrum which is roughly peaked again around the 
frequency of the fundamental mode. This was observed 
in TO] and attributed to QNM ringing. However, it is 
of completely numerical nature, with the same origin as 
the second ringing discussed in Fig. ^ when the shell 
crosses the inner boundary and leaves the grid, matter 
is artificially removed from the system, which violates 
energy conservation and produces the excitation of the 
black hole normal modes. If this happens at a relatively 
large radius such as r* = — 3M, the effect is amplified. 
Moving the inner boundary closer to the event horizon 
shifts the second QNM ringing to later times, when the 



FIG. 6: Dependence of the energy spectra on the initial lo- 
cation of the shell (k = 10), ro — 7.5M (thick solid line), 
ro — 15M (dashed line), and ro = 30M (thin solid line). 
The number of interference fringes rapidly increases with the 
initial location of the shell. 



signal is much weaker, and highly reduces the unphysi- 
cal high-frequency part of the power spectrum. There- 
fore, the energy spectrum of the GW emission that one 
could expect in a realistic astrophysical scenario, dur- 
ing anisotropic accretion of matter onto a Schwarzschild 
black hole, is most likely to be a collection of interference 
fringes covering a wide range of low frequencies, than a 
single peak at the frequency of the fundamental mode of 
the black hole. 

Next, we analyze the dependence of the energy spec- 
tra on the initial location of the shell. This is shown in 
Fig. 1^1 where we compare the spectra for three different 
initial locations, ro = 7.5M, 15M, and 30M. The num- 
ber of interference fringes rapidly increases with distance 
tq. Furthermore, the correlation between the separation 
of the fringes and the initial position of the shell is ev- 
ident: the larger the distance, the smaller the separa- 
tion between consecutive maxima. In order to explain 
these modulations in the spectrum, we can follow the 
reasoning of Lousto and Price ■ Given a GW signal 
Zo(t, r*), its Fourier transform Za[ijj,r*) is defined ac- 
cording to Eq. H23|l . For large r*, where the observer is 
located, the waveform represents only outgoing radiation, 
i.e. Zq = Zo(t — r*), which gives 

/oo 
e-*'^*Zo(< - r*)dt = e-*'^''*^(t^) , (39) 
-OO 

where 

/oo 
e-''^''Zo{u)du . (40) 
-OO 

Let us now consider another GW signal Zi whose time 
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delay with Zq is Tshift- This signal is Zi = Z{r* ,t-Tshiit), 
so that its Fourier transform reads 



Zi{uj,r* 



'Z{t-T,i,ift,r*)dt 



TABLE II: Total energy emitted in gravitational waves for 
shells accreting onto a Schwarzschild black hole from different 
distances ro. 



hiftg-^^r ^^^^ 



Thus, the Fourier transform of the signal given by the 
superposition of Zq and Zi is given by 

Z{uj,r*) = Zo{uJ,r*) + Zi{uj,r*) 

= e-*""'*A(w)(l + e-^"'^=>"f') . (41) 

Then, when computing the spectrum we have 



Z{cu, r*) ' = W\A{lu)\^ cos^ (^ujTsMit ] • (42) 



Therefore, the modulation in the frequency spectrum 
is related to the characteristic time Tghift which ac- 
counts for the delay between consecutive GW signals, 
Slu — 27r/T'shift- From the measure of the peak spacing 
du, we can thus infer the corresponding time shift. In 
the case of a radially infalling point particle, Lousto and 
Price [4^ argued that this time shift roughly coincides 
with the infalling time. They use this empirical criterion 
as a rule of thumb to predict the variation of the spectra 
as ro is increased. We have found a similar correlation 
for our extended dust shells. 

To close this section, we compute the total energy emit- 
ted in gravitational waves for shells accreting from the 
three initial locations considered previously, rg = 7.5M, 
15M, and 30M. The energy emitted is computed by 
integrating in frequency the energy spectra of Fig. El 
where the integrals are calculated using a standard trape- 
zoidal rule. As we do in the point-particle case (see Ap- 
pendix 0, it is convenient to use as a reference quantity 
the ratio 2M/fi'^E^°. Table |Il| lists the values of the en- 
ergy for the three positions ro considered. We note that 
the third value (for rg — 30M) may be slightly under- 
estimated due to some inaccuracies in the resolution of 
the power spectrum. It is worth stressing that, irrespec- 
tive of the location ro of the shell, the values reported in 
Table ^ are smaller by roughly two orders of magnitude 
than those obtained in the point-particle limit jsM li^ 
[Lousto and Price [H give {2M/fi'^)E^^ = 1.64 x lO'^ 
for ro = 30M and 1.43 x 10"^ for ro = lOM]. The reduc- 
tion we find in the total energy emitted in gravitational 
waves is a consequence of the finite size of the shells con- 
sidered in the present work. 

We note that in the numerical simulations reported 
in ^3 the estimation of the energy yielded considerably 
larger values than the ones reported here, asymptoting 
towards a third of the point-particle limit |3 as the com- 
pactness of the shell was increased. We argue that such 
a high value is overestimated, because it was affected by 
errors induced by the location of the innermost bound- 
ary of the hydrodynamics grid (r* = — 3M), as we have 
shown in Fig. \^ An indirect validation of the current 



ro/M 


(2M/^i2)£;2" 


7.5 


1.09 X IQ-"* 


15 


9.49 X IQ-^ 


30 


7.20 X IQ-^ 



estimation of energy emission comes from an inspection 
of the findings of Shapiro and Wasserman ^ . These au- 
thors compute the total energy radiated in gravitational 
waves (while here we restrict to the I = 2 multipole) 
from non-spherical dust clouds falling into a black hole 
from infinity. For any of the models considered, they 
find that the energy released in GWs is always smaller 
by at least two orders of magnitude with respect to the 
point -particle limit, with thinner clouds more efficient 
than wider ones. 



B. Neutron star simulations 

We turn now to consider the case of quadrupolar per- 
fect fluid shells accreting onto neutron stars. Two neu- 
tron star models are considered, models A and B, whose 
characteristics (mass and radius) have been described in 
Sec. Ill Al above. Initially, the quadrupolar shell is sur- 
rounded by a background fluid (an "atmosphere") sat- 
isfying the stationary and spherically symmetric Michel 
solution 40] . The initial rest mass density profile is given 
by Eq. (|37|l . where po is now the profile consistent with 
the Michel solution. As in the black hole case, the mass 
of the shell is /i = O.OIM, which corresponds to a maxi- 
mum density of pmax ~ 3.5 x 10~^ km~^ when we fix its 
width to K = 1. The (inhomogeneous) density profile of 
the atmosphere po is roughly three orders of magnitude 
lower than pmax- The shell obeys a polytropic {p = ICp^) 
EOS with K, = 0.01 km^/^ and 7 = 4/3. The initial in- 
ternal energy profile is obtained from p and p through 
the first law of thermodynamics as e = p/[(7— ^)p]- The 
shell is initially at rest at ro = 20 km, and the GW signal 
is extracted at robs = 250 km. 

We first present an overview of a typical evolution for 
both models, in order to have an immediate insight on 
the dependence of the gravitational wave signal on the 
compactness of the star. As a result of the reflecting 
boundary conditions imposed at the surface of the star 
and the existence of an extended atmosphere, the accre- 
tion process is now followed by the formation of shock 
waves which propagate off the stellar surface. As in the 
case of the black hole, the impact of the shell perturbs 
the star and triggers its quasinormal modes of pulsation. 

FigureEldisplays the time evolution of the metric func- 
tions X fc, as well as the enthalpy perturbation iJ, for 
both models. Model A (B) is presented in the left (right) 
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FIG. 7: Time evolution of the metric variables k and x ^nd of the enthalpy perturbation H for the two neutron star models 
considered. The left panels correspond to model A and the right panels to model B. The burst of gravitational radiation and 
the subsequent metric and fluid oscillations are clearly identified for both stellar models. The infalling perfect fluid shell has 
an initial compactness k = 1 and is located at a distance ro = 20 km. Note that, contrary to the previous figures, the time is 
now given in ms. 



panels. In this figure we show the all the variables and the 
logarithms of the absolute values of the metric functions, 
in order to highlight the oscillating modes of the two 
stars. Correspondingly, Fig.|Slcxhibits the time evolution 
of the Zerilli-Moncrief function for both models, in lin- 
ear and logarithmic scales, computed at every time step 
according to Eq. H18|) . In the behavior of Z, as well as in 
the time evolution of x and k shown in Fig. \7\ the three 
phases discussed in the preceding section for the black 
hole case are also visible: first, the infalling phase, when 
the bulk of the shell is evolving outside the star, grad- 
ually approaching it, which is characterized by a steady 
increase of the amplitude of the signal. This phase is 
very short as the shell is initially located very close to 



the stellar surface. Second, a burst-like peak appears in 
the GW signal, which, as found in relativistic simulations 
of gravitational core collapse , coincides with the mo- 
ment when the shell reaches the surface, creating a shock 
wave which propagates off the surface. Finally, there is 
the ringdown phase, characterized by a GW signal which 
is not exactly monochromatic as a result of the complex 
interaction between the gravitational field of the star and 
the layers of fluid captured on top of the stellar surface 
in the process of readjusting themselves to a new station- 
ary solution. The duration of the ringdown phase is now 
much longer than for the Schwarzschild black hole case 
discussed previously, as the damping time of the funda- 
mental mode of the fluid is considerably larger. We note 
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Model A 
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FIG. 8: Time evolution of the Zerilli-Moncrief function for the two stellar models considered. In close qualitative agreement 
with the results from gravitational core collapse simulations |4^ . the burst in the Zerilli-Moncrief function is associated with 
the moment of the bounce of the accreting shell at the stellar surface, and occurs at a retarded time of about 1 ms. 



that despite the waveforms obtained in our simulations 
showing a remarkable resemblance with those obtained 
in core collapse simulations (47i |. there are also impor- 
tant differences in the post-bounce phase dynamics, as 
the ringdown of the neutron star lasts for much longer 
times. In our idealized setup, with a perfect fluid, these 
pulsations are not quickly damped by the existence of a 
dense envelope surrounding the star, as happens in the 
core collapse situation. 

The bulk of the accreting matter (the center of the 
shell) reaches the stellar surface at tA ^ 0-'2 ms for model 
A and is ^ 0.14 ms for model B. At robs, where the 
observer is located, the GW signals generated by these 
events are delayed in time by 



robs -i? + 2Af log 



robs - 2Af 
R-2M 



(43) 



Therefore, the signal generated by the matter bouncing 
back at the stellar surface reaches the observer at times 
tf^ = tA+AtA « 1.06 ms andig'" = tB+Ats w 0.97 ms. 
These values are consistent with the results of Fig. [3 
where the waveforms of k and x s-re found to show bursts 
of large amplitude followed by highly damped oscillations 
sometime around these values. This observation is par- 
ticularly confirmed in the logarithmic plots of the GW 
signals. After the short-lived ringing phase lasting for 
a fraction of half a millisecond after the burst, only the 
fundamental oscillation mode of the star is visible in all 
variables plotted in Figs. Eland |H1 

The qualitative behavior found in the GW emission 
is the same for both models. Only quantitative differ- 
ences appear in the amplitudes at the maximum, which 
are systematically larger by roughly a factor of 2 for 
the more compact model (A). This difference in ampli- 



tude becomes more apparent in the energy spectra of 
the Zerilli-Moncrief function shown in Fig. O The solid 
lines in this figure are obtained by Fourier transform- 
ing the complete signal, i.e. also including the contri- 
butions from the shell infall phase which precedes the 
burst. Correspondingly, the dashed lines are obtained 
from truncated waveforms, in which we only take into 
account the contribution from the beginning of the burst 
(~ 0.9 ms) onward. We notice that, despite the short 
evolution times of our simulations, the / mode frequency 
is very well identified in the spectra, the relative differ- 
ence with respect to the values listed in Table being 
roughly of 1%. The value of the / mode frequency is 
indicated with a circle in Fig.|3 Furthermore, it is worth 
stressing that we obtain the same qualitative spectra as 
for the black hole case - a complex pattern with interfer- 
ence fringes with the addition, in the neutron star case, of 
a high peak standing at the frequency of the / mode. As 
we discussed in the preceding section for a Schwarzschild 
black hole, the comparison between the spectra obtained 
from the total and the truncated Zerilli signal shows that 
the interference fringes are produced by the interaction 
of the infalling fluid and the star. 

There are other small differences between the two stel- 
lar models. The more compact model (A) is more ef- 
ficient at high frequencies than model B. This is di- 
rectly correlated with the large amplitudes attained by 
the peaks of the Zerilli-Moncrief function in the time do- 
main (cf. Fig.|SJ). The broadband spectrum of model A 
presents one broad peak with a maximum at about 7 
kHz, which is however too low to be identified with the 
first w mode (see Tabled). As mentioned before, its ori- 
gin should be related to the motion of the fluid shell and 
its interference with the gravitational field of the star, 
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FIG. 9: Energy spectra corresponding to the time evolution of the Zerilli-Moncrief function depicted in Fig. |H| Model A (B) 
is shown on the left (right) panel. The solid lines show the spectra for the entire signal while the dashed lines are obtained 
removing that part of the waveform which corresponds to the infall phase of the shell. The spectra of the entire signal show 
characteristic interference fringes, qualitatively similar to those found in the black hole case, with a large amplitude peak 
standing at the frequency of the / mode. The value of the frequency of the / mode from Table U is indicated by a circle. The 
spectra of the truncated signals show no evidence of the interference patterns. 



or, in other words, on the reflection of the GWs pulse 
associated with the shell distribution with the "exter- 
nal" Zerilli potential. Thus, this broad feature depends 
on details of the accretion dynamics rather than on the 
intrinsic characteristics of the star. 

We note that in both spectra there appears a small 
amplitude peak at frequencies lower than that of the / 
mode for each model. This second peak is associated 
with oscillations of that part of the external fluid that 
has been gravitationally captured by the central neutron 
star as a result of accretion. In order to illustrate this 
affirmation, we plot in Fig. 1101 the energy spectra corre- 
sponding to two fluid shells which only differ on the mass 
= O.OIM and O.OOIM), falhng onto steUar model B 
from the same distance, tq = 20 km. The / mode is 
also properly excited for the less massive shell. Further- 
more, the interference fringes in the high frequency part 
of the spectra coincide to high precision for the two shells 
considered, after normalizing to the corresponding shell 
masses. However, the structure of the low frequency peak 
is the only feature of the spectra which is influenced by 
the mass of the shell. We argue that the existence of this 
unphysical low frequency peak is an artifact produced by 
the boundary conditions. Notice, once more time, that 
we model the surface of the star as a hard surface, with 
reflecting boundary conditions. In a realistic scenario, 
the accreted matter would not simply bounce at the stel- 
lar surface, but it would rather interact with the neutron 
star envelope, resulting in heating and suffering nuclear 
reactions, until it is reabsorbed by the star. 

As we did for the black hole case, let us now close this 



section by computing the total energy emitted in gravi- 
tational waves for the two stellar models considered. For 
doing so we integrate in frequency the spectra appearing 
in Fig. 1^1 The result of the integration yields 

Ef ~ 3.02 X 10"* Mqc^ , (44) 
Ef ~ 8.36 X 10"^ Mqc^ . (45) 

Model A is hence slightly more efflcient than model B 
concerning GW emission. We note in passing that these 
values are as small as those found in core collapse simula- 
tions although such a comparison only makes sense 
in qualitative terms. Finally, the mass of the shell radi- 
ated in gravitational waves is 

Ef ~ 2.15 X 10"*^ n , (46) 
Ef ~ 5.97 X lO"'^ n . (47) 



V. CONCLUSIONS 

In this paper we have presented a detailed analysis of 
the gravitational radiation induced by the anisotropic ac- 
cretion of quadrupolar (dust or perfect fluid) shells onto 
non-rotating black holes and neutron stars. The nu- 
merical framework for our simulations is based upon a 
hybrid procedure in which the linearized equations de- 
scribing metric and fluid perturbations are coupled to a 
fully nonlinear hydrodynamics code that calculates the 
motion of the accreting matter. These equations are 
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FIG. 10: Energy spectra for two fluid shells differing on the 
initial mass /i and infalling from a distance of ro = 20 km 
onto stellar model B. The solid line corresponds to /i = O.OIM 
and the dashed line to fi = O.OOlAf. The / mode excitation 
and the high frequency part of the spectra coincide to high 
precision for the two shells considered. Only the low frequency 
peak is affected by the different mass of the shell. 



integrated numerically in axisymmetry using advanced 
computational techniques. The numerical schemes devel- 
oped have proved to be stable and highly accurate. Re- 
garding the perturbation equations, the two main tech- 
nical changes with respect to previous works reported in 
the literature are the use of the Zerilli-Moncrief func- 
tion and the in-built conservation of the Hamiltonian 
constraint. The hydrodynamics equations are solved us- 
ing high-resolution shock-capturing schemes based upon 
approximate Riemann solvers. We have shown that a 
perturbative approach can be used as a very effective 
tool to understand the basic gravitational physics oper- 
ating in interesting astrophysical situations, extending 
the information which can be gained from the study of 
point-like particles infalling onto black holes or orbiting 
around them. In this context, our hybrid approach can 
be extremely useful to understand the gravitational ra- 
diation from astrophysical systems, complementing the 
whole machinery of full numerical relativity. 

The simple fluid configurations analyzed in this paper 
ascertain, however, that the effects of the extended struc- 
ture of the accreting matter are indeed highly relevant for 
gravitational wave emission. In the black hole case, most 
of the energy is released at frequencies lower than that 
of the fundamental QNM of the black hole, the spectrum 
consisting of a complex pattern, mostly produced during 
the accretion process rather than in the ringdown phase. 
Therefore, the GW emission that could be expected from 
a realistic astrophysical scenario in which accretion op- 
erates (e.g. after gravitational collapse of a massive star) 
is more likely to be a sort of interference fringe, cover- 



ing low frequencies, than a single monochromatic peak at 
the frequency of the fundamental mode of the black hole. 
This result, which went unnoticed in earlier hydrodynam- 
ical simulations by ^3 ; was already observed by Lousto 
and Price ,42] in the case of a point-like particle falling 
onto a black hole (see Appendix a situation amenable 
to semi-analytic investigation. However, we have shown 
that the appearance of interference fringes is very much 
amplified when the accreting fluid is an extended shell of 
finite size which, in turn, reduces the amount of energy 
which is released in gravitational waves to some two or- 
ders of magnitude below the point-particle value. It is 
interesting to notice that ground-based interferometric 
detectors attain the maximum sensitivity at frequencies 
considerably lower than the QNMs of stellar mass black 
holes. For this reason they are usually not considered 
as optimal sources for detection. However, to the light 
of our findings, the process of accretion appears to be 
quite more effective at frequencies about a factor of 2-3 
lower than that of the black hole fundamental mode, and 
therefore the chances of detecting gravitational wave sig- 
nals from such scenarios may be larger than expected. It 
must be stressed that the interference pattern does not 
happen only because of the GW content of the initial 
data, but it is a distinctive feature due to the extended 
size of the object, resulting from the interaction between 
the infalling matter and the backscattered waves. 

In the neutron star case the qualitative results are sim- 
ilar, but a considerable part of the energy is emitted at 
the frequency of the fundamental mode. We have shown 
that the /-mode of the star is the only one excited at 
significant levels, and that the high frequency spectrum 
is quite sensitive to the spatial distribution of the ac- 
creting matter, making the contribution of the space- 
time w-modes of the star difficult to be identified. The 
waveforms obtained show a remarkable resemblance with 
those obtained in core collapse simulations despite 
the fact that we are considering a very different scenario. 
The main difference is that in our case the ringdown of 
the neutron star lasts for much longer times: the pulsa- 
tions are not quickly damped by the existence of a dense 
envelope surrounding the star, as happens in the core 
collapse situation. 

The results reported in this paper can partly be con- 
sidered as a necessary assessment of our numerical ap- 
proach in anticipation of the study of more interesting 
astrophysical scenarios, namely the excitation of QNMs 
from perfect fluid thick accretion tori orbiting around 
compact objects (see e.g. ^HE3 references therein), 
which will be presented elsewhere ■ 
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APPENDIX A: GENERAL SOURCE TERM FOR 
THE ZERILLI-MONCRIEF EQUATION 



;> = being on a static background. Equation IIA6p be- 
comes then 

~4> + 0" + iycb' = A (-X + x") + B(-k + fc") 

+ {vC + 2C") k" - 2vCk 

+ c(^-k + k")' + L{k,x,k',x') ■ (A9) 

Next, we can write the evolution equations for k and x 
in a Schwarzschild spacetime and the Hamiltonian con- 
straint as 



In this appendix we derive the source term given in 
Eq. H21(l for a general stress-energy tensor. We use the 
normalization of the Zerilli-Moncrief function (j) given in 
Ref. |23|, although the function we evolve in the numer- 
ical simulations is rescaled a.s Z = 2(j)/X. The inhomoge- 
neous equation written using the frame derivatives nota- 
tion of Ha reads 



(Al) 



where is the potential and the source term. As 
shown by Moncrief [3l| . can be written in terms of k 
and X &s 



= A(r)x + B(r)fc + C(r)fc' , 



where 



A = 
B = 



2r^e 



2„-2b 



(A-2)r-h6M ' 

r{r\ + 2M) 
(A - 2)r + 6M ' 
C = -rAe^' . 



(A2) 

(A3) 

(A4) 
(A5) 



The derivation of the source term follows from the 
knowledge of the source terms in the evolution equations 
for X ^iid k that can be written in vacuum 27]. Using 
Eq. IIA2P in Eq. IjAip we get 



^0 + 0" + vcj)' 



= A{~x + x")-Bk 
+ [B + 2C' + vC] k" 
+ C[-(fc'r+fc"'] 
+ L(fc,x,fc',x') , 



(A6) 



where L is a linear operator acting on k, Xi and x'j 
whose explicit form is not relevant for the computation 
of the source. In fact, in Ref. [23| it is shown that all 
terms which are linear in the fields x ^-^d k and their 
first-order spatial frame derivatives merge together to 
build the Zerilli potential. The frame derivatives do not 
commute Hal, 



so that 



{ky-{k'j=-uf 



{k'y^ (k)' + 2vk 



(A7) 



(A8) 



-X + X" ^L{k,x.x'} + T^ (AlO) 
-k + k" {x, k} - {2W + v) k' + Tk (All) 
k" ^L{x,k,x',k'} + Tn, (A12) 

with r^, Tfe and T-^ being the sources induced by the 
matter distribution. From these equations one obtains 



-k + k' 



^L{k,x,k',x'}-i2W + ,^)Tn + U , 

(A13) 



so that the source of the Zerilli equation is found to be 
= AT^ + {B + 2vC)Tk 

+ [2C' -2{W + iy)C]Tn + CTl: ■ (A14) 
In polar radial coordinates the sources are, explicitly. 



T-H = -8ttS' 



Ti. , 



(A15) 
(A16) 



Tfc = 87r<i - 



e 111 



-2b 



r r 
A + 2 66-26 



^2 



(A17) 



where Sn and have been defined in Sec. Ill ('I When 
Eqs. HA15|) - (jA17p and the definitions of the coefficients 
(|J3|) - (|X5|l are replaced in Eq. HA14|) . we get = A5^/2, 
with Sz given by Eq. |(5lJ . 

The projections of the source stress-energy tensor onto 
the basis of the spherical harmonics is accomplished as 
follows. From the orthogonality properties of the har- 
monics. 



dnY;„^^Yi,,r^',br''^xSu'S„ 



I 



dnz. 



I'm' r,ab _ A (A - 2) 



ab 



(A18) 
(A19) 
(A20) 
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FIG. 11: Gravitational waveforms and energy spectra for a point particle falling radially onto a Schwarzschild black hole from 
different distances. The left panels show the waveforms and the right panel the corresponding spectral energy distributions. 
Excellent agreement is found with the results of Lousto and Price |4^. See text for further details. 



we obtain the coefficients of the expansion of tuii, 
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where, following Regge and Wheeler |32j, we have defined 
the functions 



^Im = 2im [coti?y/„ - YLxA 



W, 



(A26) 
(A27) 



In the general case where t^i, describes a complex source 
corresponding to a general distribution of matter evolv- 
ing dynamically, these integrals have to be evaluated nu- 
merically. The analytic computation can be done only 
in some simple cases, such as the particular case when 
the source is a test-mass body moving along a geodesic 
of Schwarzschild spacetime 



APPENDIX B: POINT-LIKE PARTICLES 
RADIALLY FALLING ONTO BLACK HOLES. 

In this appendix we reexamine the simplified scenario 
of a point particle radially falling onto a Schwarzschild 
black hole. This is done with two purposes: ffi'st, to 
test our perturbative numerical code with previous works 
and, second, to analyze the similarities and differences 
with the case of accretion of extended fluid shells onto 
black holes. The emission generated by infalling parti- 
cles has been extensively studied in the past. The seminal 
calculation of the GW emission when a test particle falls 
from infinity Q was later extended to non-radial trajec- 
tories by Detweiler and Szedenits [S^l- In both cases the 
analysis was done in the frequency domain. A frequency 
component treatment based on Laplace transforms was 
also employed by Lousto and Price study of 

the emission from particles falling from finite distances. 
On the other hand, a treatment of the same problem in 
the time domain hasjust recently been approached by 
Mart el and Poisson |3a|. 

In order to test our numerical code we should be able to 
reproduce those results reported in Refs. | ,38i . i43j within 
the current time domain approach. A way to deal with 
a (5-like source and with a discontinuous Zerilli-Moncrief 
function at the location of the particle was developed by 
Lousto and Price 0|> ^^'^ later successfully applied in 
Ref. [3^ following a time domain approach. However, we 
found it convenient and accurate enough for our purposes 
to represent the S function in the particle source terms 
by a narrow Gaussian written as 
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FIG. 12: Effects of tfie gravitational wave content of the ini- 
tial data on the energy spectra of a point-Hke particle falling 
radially onto a black hole from Ro = 15M. By eliminating 
the initial contribution of GWs the amplitude of the bumps is 
strongly reduced. Qualitatively, however, the effect is always 
present. 



with CT ^ 1. For a given resolution, we fix cr = 2Ar*, 
so that the smaller Ar* , the better the approximation of 
the 6 function is. The Zerilli-Moncrief equation is solved 
on an evenly spaced grid using a standard three-level 
leapfrog scheme. We consider particles initially at rest 
falling from a finite distance Rq to compare with the 
results of Ref. Since the particle is falling along 

the z axis, the system is axisymmetric and the only non- 
vanishing contribution for any I is the m — one. The 
source terms are specified accordingly and we consider 
explicitly just I — 2. 

As a test of our numerical method, we consider the 
particle falling from Rq = lOM, 15M, and 20M, with 



the same initial setup of Ref. [43 (i-e. including some 
initial GW content) and compute the evolution of the 
Zerilli-Moncrief function and its corresponding energy 
spectrum. Figure shows the results of these simula- 
tions. The left panels show the temporal evolution of Z 
normalized to the particle mass, while the right panel ex- 
hibits the corresponding energy spectra. The waveforms 
show a good agreement in amplitude and shape with 
those of Ref. [43 ■ The energy spectra for Rq = lOM, 
15M, and 20M must be compared with Figs. 6(b), 4(b), 
and 4(c) of Ref. respectively. The spectra show the 
same localized bumps due to the interference between the 
initial GW pulse and the GWs emitted by the particle 
during its motion. We notice quite good agreement for 
Rq = lOM and Rq = 15M, while some small differences 
are found for Rq = 20M. 

Next, we analyze the contribution of the initial data on 
the power spectra. We have studied how the spectrum 
changes when the initial GW contribution is eliminated 
from the evolution. This is accomplished in the same 
way we used for the extended shells; that is, the particle 
is frozen at its location until the initial pulse has gone 
from the numerical domain, after which the evolution 
starts. Figure 1121 compares the energy spectra emitted 
by a particle falling from Rq = 15M with (dashed line) 
and without (solid line) the initial GW contribution. The 
signal is extracted at r = 500M in the two cases. As sug- 
gested in a more general scenario in Ref. js^ , our results 
confirm that in the test particle case the bumps in the 
spectrum are mainly due to the spurious contribution of 
the radiation in the initial data. In fact, when removing 
the initial GW content, the amplitude of the bumps is 
reduced. Some modulation is, nevertheless, still present 
and the spectrum does not fully correspond to that of 
a pure QNM ringdown signal. Martel and Poisson [s^ 
argued that this is to be interpreted as an interference 
effect as well, but between the waves emitted by the par- 
ticle during its motion and those previously emitted and 
backscattered by the potential. 
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